A Pressure-Based Fully-Coupled Flow Algorithm for the Control Volume Finite Element Method

Abstract

A pressure-based fully-coupled flow algorithm is developed for the control volume finite element method (CVFEM), which covers speeds up to the transonic regime. The CVFEM is used because it presents a number of advantages as compared to the popular cell-centered finite volume method (CCFVM), while retaining the properties of the finite volume method (FVM) in terms of flux conservation and numerical integration simplicity. The implementation presents a novel weak formulation of Dirichlet boundary conditions to resolve the disadvantages emerging from the strong formulation, by mimicking the methods followed in the CCFVM. Derivation and implementation details are presented, and a number of test cases are used to evaluate the accuracy and performance of this approach.

1. Introduction

Computational fluid dynamics (CFD) is a critical tool for a wide range of industries. As such, it continues to be the focus of continuous research and development, with the aim of increasing robustness, performance, and general capabilities. To this date, the majority of CFD solvers have been based on the CCFVM as it is inherently conservative and presents a simple numerical framework. It also allows for a relatively easy implementation applied directly over a mesh with no additional layers of geometric manipulations. However, while the FVM can handle all types of meshes, its accuracy and convergence characteristics are greatly affected by the mesh quality, as defined in terms of high aspect ratio and severe non-orthogonality, in addition to the topology of mesh elements [1,2]. The CVFEM, known also as vertex-centered FVM approach [3,4], can handle meshes with suboptimal properties with less difficulty. This is mainly due to its use of shape functions that yield improved accuracy for the computed gradients, and support a much more comprehensive linearization of field fluxes. This in turn allows for a more implicit discretization, and thus a more robust solution process. The CVFEM is also fully conservative and that is equally suited to the simulation of fluid flow problems. It is, however, computationally more expensive than the CCFVM in terms of memory and programming complexity. This has played a role in the popularity of the CCFVM methods, even though they have both emerged at roughly the same time, but, nevertheless, and in view of the above properties, the last decade has witnessed a renewed interest in the CVFEM; the method has been the subject of advancing efforts in three-dimensional advection–diffusion problems, incompressible and compressible flows, two-phase flows, flow in porous media, and also in adapting it for high-order schemes [5,6,7,8,9,10,11,12,13,14]. It is also worth mentioning the extensive use of this method in disciplines such as nanofluids, nanoparticles, and nanomaterials [15,16,17], and realization in heat transfer applications, pipeline design, and reservoir engineering [18,19,20,21,22,23].

Furthermore, there are now a number of open-source libraries [24] that provide the initial building blocks for the implementation of advanced solution algorithms and numerical techniques for the CVFEM. It is interesting, however, to mention an additional advantage of the CVFEM approach in that it is not bound to be second-order accurate, but can be adapted for higher-order accuracy, regardless of the mesh refinement level. This can be extremely useful for low-quality meshes, and it may also be the foundation of a new generation of CFD solvers, because it provides the means of reaching high-order accuracy with minimal cost of meshing time. Nonetheless, this advantage is beyond the outcomes of the current work, which is intended to be second-order accurate, and is to be very possibly part of an upcoming development effort.

In this work, the fully-coupled algorithm [25,26,27,28,29,30,31] is reformulated for the CVFEM. The algorithm is implemented as a basis for an in-house CFD solver, aiming at solving complex engineering fluid flow problems ranging from incompressible up to transonic-compressible regimes. The fully-coupled algorithm is crucial for achieving good robustness and performance, and is more critical in this case since an iteration of the CVFEM is more expensive than a similar iteration with the CCFVM. The challenges to address are numerous, starting with the treatment of the boundary conditions that had to be completely rethought, since boundary nodes are now part of the computational domain and are involved in the solution vector [32,33], but also in determining how the velocity–pressure coupling can be expressed in terms of shape functions. In what follows, the governing equations are introduced, and then the CVFEM method is briefly reviewed. The implementation of the velocity–pressure coupling is then introduced in detail, followed by a derivation of a number of boundary conditions. The method is then tested on a number of flow problems for incompressible and compressible flows, and the performance, robustness, and accuracy of the method are evaluated.

2. Governing Equations

The flow equations are mainly described by the mass, momentum, and energy conservation equations. The energy equation is presented in its specific total enthalpy form. The equations are, respectively, as follows:

$$\begin{array}{cc}\hfill \frac{\partial \rho }{\partial t}+\nabla ·\left(\rho \mathbf{v}\right)& =0\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill \frac{\partial \rho \mathbf{v}}{\partial t}+\nabla ·\left(\rho \mathbf{v}\mathbf{v}\right)& =\nabla ·\tau -\nabla p\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill \frac{\partial \rho h_0}{\partial t}+\nabla ·\left(\rho \mathbf{v}{h}_0\right)& =\nabla ·({\lambda }_{eff}\nabla T)+Q^{h_0}\hfill \end{array}$$

(3)

where $\tau$ in Equation (2) is the stress tensor which includes laminar and turbulent stresses. Assuming a Boussinesq assumption, $\tau$ is expressed as

$$\tau ={\mu }_{eff}\left(\nabla \mathbf{v}+\nabla {\mathbf{v}}^T\right)-\frac{2}{3}{\mu }_{eff}(\nabla ·\mathbf{v})\mathbf{I}-\frac{2}{3}\rho k\mathbf{I}$$

(4)

${\mu }_{eff}$ and ${\lambda }_{eff}$ are effective viscosity and effective thermal conductivity, respectively, defined as

$${\mu }_{eff}=\mu +{\mu }_t,{\lambda }_{eff}=\lambda +{\lambda }_t$$

$Q^{h_0}$ in the total enthalpy equation consists of energy sources such as viscous dissipation. An auxiliary equation, i.e., the ideal gas law, is required to relate density to pressure and temperature, stated here:

$$\rho =\frac{p}{RT}$$

(5)

The system of equations above is closed by the SST $k-\omega$ turbulence model [34]:

$$\begin{array}{cc}\hfill \frac{\partial \rho k}{\partial t}+\nabla .\left(\rho \mathbf{v}k\right)& ={\overline{\tau }}_{ij}\nabla \mathbf{v}-{\beta }^{*}\rho \omega k+\nabla .\left[(\mu +{\sigma }_k{\mu }_t)\nabla k\right]\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill \frac{\partial \rho \omega }{\partial t}+\nabla .\left(\rho \mathbf{v}\omega \right)& =\frac{\rho \gamma {\overline{\tau }}_{ij}\nabla \mathbf{v}}{{\mu }_t}-\beta \rho {\omega }^2+\nabla .\left[(\mu +{\sigma }_{\omega }{\mu }_t)\nabla \omega \right]+2(1-F_1)\frac{\rho {\sigma }_{\omega 2}}{\omega }\nabla k\nabla \omega \hfill \end{array}$$

(7)

The model’s parameters, blending functions, and variables are retrieved from Ref. [34], and may not be stated here.

3. The Control Volume Finite Element Method

CVFEM is a combination of the FVM and the finite element method (FEM). The standard FVM utilizes a cell-centered scheme, with control volume and the element set as identical, and the variables stored at the control volume centroid. This is shown in Figure 1a with control volume shaded in gray and the variable set at the element centroid represented by a red point. The CVFEM is a vertex-centered scheme that stores the variable of interest at the vertices (nodes) of the mesh elements. The control volume is based on the dual mesh and is defined around the element nodes. The control volumes are thus not identical with the mesh elements; this is portrayed in Figure 1b, where the control volume, gray-shaded, is centered at the vertex. In the CVFEM, control volumes are also constructed at boundary nodes, because variables are stored on nodes; hence, these nodes are part of the solution vector. Boundary control volumes for CVFEM, on the other hand, are shown in Figure 2, where it is obvious that control volumes only partially encapsulate the nodes.

Integration points ($ip$) are defined for the CVFEM on the surface of the dual volume, just as for the CCFVM, where they are also defined at the element face centroid. Integration points, shown at the centroids of the sub-control surfaces in Figure 3a, are used for the numerical integration of surface integrals with the use of arbitrary interpolation functions, unique to each element topology. The value of a generic quantity ${\varphi }_{ip}$ at an integration point is therefore computed based on the interpolation weights corresponding to the nodal values of ${\varphi }_k$.

$${\varphi }_{ip}=\sum _{k=1}^n{\left.N_k\right|}_{ip}{\varphi }_k$$

(8)

n is the number of nodes of the element containing $ip$ and ${\varphi }_k$ the value at the node of index k straddling $ip$, whereas the calculation of the gradient of the quantity $\varphi$ at the integration point is given by

$$\nabla {\varphi }_{ip}=\sum _{k=1}^n{\left.\nabla N_k\right|}_{ip}{\varphi }_k$$

(9)

While interior integration points involve all nodes of the element, only the nodes of a boundary face are required to interpolate the nodal values to the boundary integration point. This is shown in Figure 3a. That is due to the definition of the shape function to be zero at edge-opposite nodes if the integration point is on an edge. Therefore, the value of $\varphi$ at a boundary integration point is

$${\varphi }_{bip}=\sum _{k=1}^{n^s}{\left.N_k\right|}_{bip}{\varphi }_k$$

(10)

where $n^s$ is the number of boundary (side) nodes straddling $bip$. However, this is not the case when computing a gradient at a boundary integration point, because there do exist non-zero values for the derivative of shape functions for interior nodes contributing to boundary integration points (Figure 3c). Therefore, the gradient of $\varphi$ at $bip$ is computed as

$$\nabla {\varphi }_{bip}=\sum _{k=1}^n{\left.\nabla N_k\right|}_{bip}{\varphi }_k$$

(11)

4. The Pressure-Based Coupled Solution Procedure

In a fully-coupled velocity–pressure algorithm, the momentum and mass pressure equations are solved simultaneously. The structure of the system is as follows:

$$\begin{array}{ccc}\mathrm{Momentum}\mathrm{Equation}:& & \\ \mathrm{Mass}\mathrm{Equation}:& \end{array}\left[\begin{array}{cc}{\mathbf{a}}_C^{\mathbf{vv}}& {\mathbf{a}}_C^{\mathbf{v}p}\\ {\mathbf{a}}_C^{p\mathbf{v}}& a_C^{pp}\end{array}\right]\left[\begin{array}{c}{\mathbf{v}}_i\\ p_i\end{array}\right]+\sum _{F,j}\left[\begin{array}{cc}{\mathbf{a}}_F^{\mathbf{vv}}& {\mathbf{a}}_F^{\mathbf{v}p}\\ {\mathbf{a}}_F^{p\mathbf{v}}& a_F^{pp}\end{array}\right]\left[\begin{array}{c}{\mathbf{v}}_j\\ p_j\end{array}\right]=\left[\begin{array}{c}{\mathbf{b}}_C^{\mathbf{v}}\\ b_C^p\end{array}\right]$$

(12)

where the subscript C appearing in Equation (12) denotes a control volume associated with the node of index i, and F is any neighboring control volume to C, associated with a node of index j. There are, however, two categories of neighboring control volumes F; the first category includes those that share an integration point with C, while the other category includes those that do not share an integration point with C (Figure 4). The only tangible difference which arises between them is related to the advection fluxes. Neighboring control volumes that do not share an integration point with C are not involved in the discretization of advection fluxes. This will be highlighted in the following section. Superscripts * and ∘ will denote the latest available value and the previous time step value, respectively, while the notation $\left|\right|a,b\left|\right|$ will be used to indicate a maximum between two quantities a and b.

4.1. Discretized Momentum Conservation Equation

The discrete momentum conservation equation appearing in the coupled system (Equation (12)) is restated here:

$${\mathbf{a}}_C^{\mathbf{vv}}{\mathbf{v}}_i+\sum _{F,j}{\mathbf{a}}_F^{\mathbf{vv}}{\mathbf{v}}_j+{\mathbf{a}}_C^{\mathbf{v}p}{p}_i+\sum _{F,j}{\mathbf{a}}_F^{\mathbf{v}p}{p}_j={\mathbf{b}}_C^{\mathbf{v}}$$

(13)

The coefficients in Equation (13) are obtained by applying the following discretization method for the various terms: the implicit backward Euler for transient, the upwind scheme for the advection, a symmetric interpolation for diffusion, and Gauss integration of a control volume for pressure gradient. For the diffusion and pressure gradient, the interpolation to the integration points is performed using the shape functions. This yields the following expressions:

$$\begin{array}{cc}\hfill {\mathbf{a}}_C^{\mathbf{vv}}& =\left\{\underset{\mathrm{transient}}{\underbrace{\frac{{\rho }_i^{*}{V}_C}{\Delta t}}}+\sum _{ip}\left(\underset{\mathrm{stress}-\mathrm{part}1}{\underbrace{-{\mu }_{ip}{\left.\nabla N_i\right|}_{ip}·{\mathbf{S}}_{ip}}}+\underset{\mathrm{advection}}{\underbrace{\left|\right|{\dot{m}}_{ip}^{*},0\left|\right|}}\right)\right\}\mathbf{I}+\sum _{ip}\underset{\mathrm{stress}-\mathrm{part}2}{\underbrace{\left(-\frac{1}{3}{\mu }_{ip}{\left.\nabla N_i\right|}_{ip}{\mathbf{S}}_{ip}\right)}}\hfill \\ \hfill {\mathbf{a}}_F^{\mathbf{vv}}& =\left(\underset{\mathrm{stress}-\mathrm{part}1}{\underbrace{-{\mu }_{ip}{\left.\nabla N_j\right|}_{ip}·{\mathbf{S}}_{ip}}}-\underset{\begin{array}{c}\mathrm{advection}\\ 0\mathrm{if}ip\mathrm{not}\mathrm{shared}\mathrm{by}F\end{array}}{\underbrace{\left|\right|-{\dot{m}}_{ip}^{*},0\left|\right|}}\right)\mathbf{I}+\underset{\mathrm{stress}-\mathrm{part}2}{\underbrace{\left(-\frac{1}{3}{\mu }_{ip}{\left.\nabla N_j\right|}_{ip}{\mathbf{S}}_{ip}\right)}}\hfill \\ \hfill {\mathbf{a}}_C^{\mathbf{v}p}& =\sum _{ip}\underset{\mathrm{pressure\; gradient}}{\underbrace{{\left.N_i\right|}_{ip}{\mathbf{S}}_{ip}}}\hfill \\ \hfill {\mathbf{a}}_F^{\mathbf{v}p}& =\underset{\mathrm{pressure\; gradient}}{\underbrace{{\left.N_j\right|}_{ip}{\mathbf{S}}_{ip}}}\hfill \\ \hfill {\mathbf{b}}_C^{\mathbf{v}}& =\underset{\mathrm{transient}}{\underbrace{\frac{{\rho }_i^{\circ }{V}_C}{\Delta t}{\mathbf{v}}_i^{\circ }}}-\sum _{ip}\underset{\mathrm{stress}-\mathrm{part}3}{\underbrace{\frac{2}{3}{\rho }_{ip}^{*}{k}_{ip}^{*}{\mathbf{S}}_{ip}}}\hfill \end{array}$$

(14)

4.2. Discretized Mass Conservation Equation

The pressure equation appearing in Equation (12) can be written as

$${\mathbf{a}}_C^{p\mathbf{v}}{\mathbf{v}}_i+\sum _{F,j}{\mathbf{a}}_F^{p\mathbf{v}}{\mathbf{v}}_j+a_C^{pp}{p}_i+\sum _{F,j}{a}_F^{pp}{p}_j=b_C^p$$

(15)

and the coefficients are defined as follows:

$$\begin{array}{cc}\hfill {\mathbf{a}}_C^{p\mathbf{v}}& =\sum _{ip}\underset{\mathrm{mass\; divergence}}{\underbrace{{\rho }_{ip}^{*}{\left.N_i\right|}_{ip}{\mathbf{S}}_{ip}}}\hfill \\ \hfill {\mathbf{a}}_F^{p\mathbf{v}}& =\underset{\mathrm{mass\; divergence}}{\underbrace{{\rho }_{ip}^{*}{\left.N_j\right|}_{ip}{\mathbf{S}}_{ip}}}\hfill \\ \hfill a_C^{pp}& =\underset{\mathrm{transient}}{\underbrace{\frac{{\left(\frac{\partial \rho }{\partial p}\right)}_{p_i^{*}}{V}_C}{\Delta t}}}+\sum _{ip}\underset{\mathrm{advection}-\mathrm{like}}{\underbrace{\left|\right|\frac{{\dot{m}}_{ip}^{*}}{{\rho }_{ip}^{*}},0\left|\right|{\left.N_i\right|}_{ip}{\left(\frac{\partial \rho }{\partial p}\right)}_{p_i^{*}}}}-\sum _{ip}\underset{\mathrm{diffusion}-\mathrm{like}}{\underbrace{{\rho }_{ip}^{*}\overline{{\mathbf{D}}_{ip}^{\mathbf{v}}}{\left.\nabla N_i\right|}_{ip}·{\mathbf{S}}_{ip}}}\hfill \\ \hfill a_F^{pp}& =-\underset{\begin{array}{c}\mathrm{advection}-\mathrm{like}\\ 0\mathrm{if}ip\mathrm{not}\mathrm{shared}\mathrm{by}F\end{array}}{\underbrace{\left|\right|-\frac{{\dot{m}}_{ip}^{*}}{{\rho }_{ip}^{*}},0\left|\right|{\left.N_j\right|}_{ip}{\left(\frac{\partial \rho }{\partial p}\right)}_{p_j^{*}}}}-\underset{\mathrm{diffusion}-\mathrm{like}}{\underbrace{{\rho }_{ip}^{*}\overline{{\mathbf{D}}_{ip}^{\mathbf{v}}}{\left.\nabla N_j\right|}_{ip}·{\mathbf{S}}_{ip}}}\hfill \\ \hfill b_C^p& =\underset{\mathrm{transient}}{\underbrace{\frac{{\rho }_i^{\circ }{V}_C}{\Delta t}}}\hfill \end{array}$$

(16)

For an ideal gas, the partial derivative $\frac{\partial \rho }{\partial p}$ is given by

$$\frac{\partial \rho }{\partial p}=\frac{1}{RT}$$

(17)

$\overline{{\mathbf{D}}_{ip}^{\mathbf{v}}}$ is an adaptive tensor interpolated to an integration point from the straddling nodal values:

$$\overline{{\mathbf{D}}_{ip}^{\mathbf{v}}}=\sum _{k=1}^n{\left.N_k\right|}_{ip}{\mathbf{D}}_k^{\mathbf{v}}$$

(18)

and for a node of index i, ${\mathbf{D}}_i^{\mathbf{v}}$ can be expressed such as:

$${\mathbf{D}}_i^{\mathbf{v}}=\left[\begin{array}{ccc}{D}_i^{v^x}& 0& 0\\ 0& D_i^{v^y}& 0\\ 0& 0& D_i^{v^z}\end{array}\right]=\left[\begin{array}{ccc}\frac{V_C}{a_C^{v^x{v}^x}}& 0& 0\\ 0& \frac{V_C}{a_C^{v^y{v}^y}}& 0\\ 0& 0& \frac{V_C}{a_C^{v^z{v}^z}}\end{array}\right]$$

(19)

$v^x$, $v^y$ and $v^z$ are the Cartesian components of $\mathbf{v}$.

4.3. Discretized Energy Conservation Equation

The energy conservation equation is a trivial scalar transport equation with no special treatments. Considering an ideal gas, and assuming a constant thermal conductivity, the discrete form is given by

$$a_C^{h_0}{h}_{0,i}+\sum _{F,j}{a}_F^{h_0}{h}_{0,j}=b_C^{h_0}$$

(20)

where the coefficients are defined as

$$\begin{array}{cc}\hfill a_C^{h_0}& =\underset{\mathrm{transient}}{\underbrace{\frac{{\rho }_i^{*}{V}_C}{\Delta t}}}+\sum _{ip}\underset{\mathrm{advection}}{\underbrace{\left|\right|{\dot{m}}_{ip}^{*},0\left|\right|}}-\sum _{ip}\underset{\mathrm{diffusion}}{\underbrace{\frac{{\lambda }_{ip}}{c_{p,ip}}{\left.\nabla N_i\right|}_{ip}·{\mathbf{S}}_{ip}}}\hfill \\ \hfill a_F^{h_0}& =-\underset{\begin{array}{c}\mathrm{advection}\\ 0\mathrm{if}ip\mathrm{not}\mathrm{shared}\mathrm{by}F\end{array}}{\underbrace{\left|\right|-{\dot{m}}_{ip}^{*},0\left|\right|}}-\underset{\mathrm{diffusion}}{\underbrace{\frac{{\lambda }_{ip}}{c_{p,ip}}{\left.\nabla N_j\right|}_{ip}·{\mathbf{S}}_{ip}}}\hfill \\ \hfill b_C^{h_0}& =\underset{\mathrm{transient}}{\underbrace{\frac{{\rho }_i^{\circ }{V}_C}{\Delta t}{h}_{0,i}^{\circ }}}+\underset{\mathrm{source}}{\underbrace{Q^{h_0}{V}_C}}\hfill \end{array}$$

(21)

5. Boundary Conditions

As opposed to the CCFVM, where Dirichlet boundary conditions are imposed at boundary faces, they are rather imposed at boundary nodes in FEM. Being part of the solution vector, the Dirichlet values of these boundary nodes are directly set in the vector, a process known as the strong formulation of Dirichlet boundary conditions. This formulation, however, gives rise to an ambiguity at nodes shared by two different adjacent boundaries. It is very common, for instance, to have a Dirichlet boundary adjacent to a Neumann boundary (Figure 5). In this case, the Neumann boundary will have no effect on the common node, a fact which results in a non-physical behavior. Degraded robustness is also reported for a strong formulation of the Dirichlet boundary [32,33]. This same dilemma exists also for the CVFEM; however, being an FVM-based method, this allows for a weak formulation of Dirichlet boundary conditions in a very similar way to the CCFVM, where fluxes are computed at boundary faces and added as a contribution to the associated coefficients. This also implies that the boundary nodes will be influenced by interior node values, and, thus, they are not expected to have values exactly equal to the Dirichlet condition, but close enough at convergence state. In this work, a novel weak formulation of the Dirichlet boundary conditions is introduced so as to replace the strong formulation, aiming at solving the aforementioned ambiguity as well as overcoming the degraded robustness problem. In what follows, a brief description of a number of boundary conditions is stated in a weak formulation, that is, similar to the way they are incorporated in a typical the CCFVM.

5.1. Boundary Conditions for the Momentum Conservation Equation

5.1.1. No-Slip Wall

With a zero mass flux, there is no contribution from the advection term to the coefficients of a wall. However, the shear stress term contributes to the coefficients in the wall’s tangential direction. Hence, the shear force at the wall boundary integration point ($bip$) can be written as

$${\mathbf{F}}_{bip}={\left({\mathbf{F}}_{bip}\right)}_{\Vert }={\left({\tau }_{bip}\right)}_{\Vert }{S}_{bip}$$

(22)

where ${\left({\tau }_{bip}\right)}_{\Vert }$ is the tangential stress at $bip$ and is given by

$${\left({\tau }_{bip}\right)}_{\Vert }=-{\mu }_{bip}{\left.\frac{\partial {\mathbf{v}}_{\Vert }}{\partial d_{\perp }}\right|}_{bip}$$

(23)

${\mathbf{v}}_{\Vert }$ is the tangential velocity, and $d_{\perp }$ is the normal distance from the wall. The velocity parallel to the wall ${\mathbf{v}}_{\Vert }$ is computed by a vector projection with the unit normal vector of the wall $\mathbf{n}$:

$${\mathbf{v}}_{\Vert }=\mathbf{v}-(\mathbf{v}·\mathbf{n})\mathbf{n}$$

(24)

Moreover, the pressure at the wall can be interpolated from the boundary control volume node values.

5.1.2. Inlet

The simplest inlet boundary condition is a specified velocity. This completely defines the boundary mass flux:

$${\dot{m}}_{bip}={\rho }_{bip}{\left.{\mathbf{v}}_{spec}\right|}_{bip}·{\mathbf{S}}_{bip}$$

(25)

It is also worth mentioning here the implementation of an inlet boundary condition where the total pressure $p_0$ is specified, in addition to the direction of the flow. For compressible flows, $p_0$ is expressed as

$$p_0=p{\left(1+\frac{\gamma -1}{2}\frac{v^2}{\gamma RT}\right)}^{\frac{\gamma -1}{\gamma }}$$

(26)

v and p are obviously unspecified, but the method followed in this work assumes that p is implicitly specified by the formula:

$$p=p_0{\left(1+\frac{\gamma -1}{2}\frac{v^2}{\gamma RT}\right)}^{-\frac{\gamma -1}{\gamma }}$$

(27)

The flow direction is also implicitly implemented in the advection term, such that the advected velocity is given by

$${\mathbf{v}}_{bip}=\frac{\left({\sum }_{k=1}^{n^s}{\left.N_k\right|}_{bip}{\mathbf{v}}_k\right)·{\mathbf{n}}_{bip}}{{\mathbf{n}}_{bip}·{\mathbf{e}}_v}{\mathbf{e}}_v$$

(28)

where ${\mathbf{e}}_v$ is the specified flow direction. The pressure gradient uses the prescribed pressure values at the boundary, so this will only contribute to the right-hand side of the system.

5.1.3. Outlet

At an outlet, the pressure is usually prescribed. The shear stress term assumes a fully developed flow, and therefore no stress, normal to the boundary. The pressure gradient uses the prescribed pressure values at the boundary, so this will only contribute to the right-hand side of the system.

5.1.4. Opening

An opening is a boundary which allows the flow to cross in either directions. In case of outflow, the boundary is treated exactly as an outlet with specified pressure, while in case of inflow, it is treated exactly as an inlet with a specified total pressure. Hence, the modifications to the coefficients are the same. The solver will obviously check the flow direction at each boundary integration point in order to apply the suitable discretization.

5.1.5. Symmetry Plane

There is no mass flux through a symmetry plane; therefore, there is no contribution from the advection term to the coefficients. The shear stress contributes to the coefficients in the symmetry’s normal direction. Hence, the wall force at a boundary integration point ($bip$) can be written as

$${\mathbf{F}}_{bip}={\left({\mathbf{F}}_{bip}\right)}_{\perp }={\left({\tau }_{bip}\right)}_{\perp }{S}_{bip}$$

(29)

where ${\left({\tau }_{bip}\right)}_{\perp }$ is the normal stress and can be approximated, according to Ref. [35], by

$${\left({\tau }_{bip}\right)}_{\perp }\approx -2{\mu }_{bip}{\left.\frac{\partial {\mathbf{v}}_{\perp }}{\partial d_{\perp }}\right|}_{bip}$$

(30)

${\mathbf{v}}_{\perp }$ is the normal velocity, computed as

$${\mathbf{v}}_{\perp }=(\mathbf{v}·\mathbf{n})\mathbf{n}$$

(31)

Similarly to the case of a wall, the derivative in Equation (30) can be computed by the use of shape functions’ gradients. In addition, the pressure at a symmetry plane usually assumes the average of the interior.

5.2. Boundary Conditions for the Mass Conservation Equation

For incompressible and subsonic compressible flows, two types of boundary conditions mainly exist for the pressure equation. The first type is that which occurs on a wall or an inlet where velocity is specified, and the second type is when the pressure is specified, such as in the case of an outlet. In the first type, the mass flux at any sub-control surface is known, while in the second, since the pressure is specified, the situation is considered to be a Dirichlet-type boundary condition.

5.2.1. No-Slip Walls

As mentioned earlier, there is no mass flux through a wall, so the mass flux at any wall sub-control surface is specified as ${\dot{m}}_{bip}=0$. Therefore, there is no contribution to the matrix coefficients.

5.2.2. Inlet

For a specified velocity inlet and due to the fact that weak formulations are utilized, the mass flux is not simply computed from the specified velocity, but in fact is determined from the Rhie–Chow formula, as follows:

$${\dot{m}}_{bip}={\rho }_{bip}\left[{\left.{\mathbf{v}}_{spec}\right|}_{bip}-\overline{{\mathbf{D}}_{bip}^{\mathbf{v}}}\left(\nabla p_{bip}-\overline{\nabla p_{bip}^{*}}\right)\right]·{\mathbf{S}}_{bip}$$

(32)

where $\nabla p_{bip}$ is evaluated as

$$\nabla p_{bip}=\sum _{k=1}^n{\left.\nabla N_k\right|}_{bip}{p}_k$$

(33)

If the total pressure is specified at the inlet, the treatment resembles that of a specified pressure at an outlet, which will come next. The only difference is that the specified pressure is the implicitly specified pressure stated in Equation (27) previously.

5.2.3. Outlet

At the outlet the pressure is specified and the mass flow rate is computed using the Rhie–Chow formula, as follows:

$${\dot{m}}_{bip}={\rho }_{bip}\left[{\mathbf{v}}_{bip}-\overline{{\mathbf{D}}_{bip}^{\mathbf{v}}}\left(\nabla p_{bip}-\overline{\nabla p_{bip}^{*}}\right)\right]·{\mathbf{S}}_{bip}$$

(34)

where ${\mathbf{v}}_{bip}$ and $\nabla p_{bip}$ are evaluated as

$$\begin{array}{cc}\hfill {\mathbf{v}}_{bip}& =\sum _{k=1}^{n^s}{\left.N_k\right|}_{bip}{\mathbf{v}}_k\hfill \\ \hfill \nabla p_{bip}& =\sum _{k=1}^n{\left.\nabla N_k\right|}_{bip}{p}_k\hfill \end{array}$$

(35)

5.2.4. Opening

The modifications to the mass-related coefficients of the matrix are the same as those applied for a specified pressure outlet in the case of outflow, and specified total pressure inlet in the case of inflow.

5.2.5. Symmetry Plane

Similar to the wall boundary condition, the mass flux through the boundary is zero, so this leaves the matrix untouched.

5.3. Boundary Conditions for the Energy Conservation Equation

5.3.1. No-Slip Wall

No contributions to the wall boundary nodes exist, except for the case of a specified temperature or a specified heat flux. If a heat flux is set at the wall, then only the right-hand side of the system is modified.

5.3.2. Inlet

Temperature is commonly specified at an inlet. Therefore, there are contributions from advection and diffusion terms to the coefficients of the boundary nodes.

5.3.3. Outlet

For the energy conservation equation, the temperature usually has a boundary condition of zero-gradient. Therefore, there is only an advection contribution to the diagonal.

5.3.4. Opening

The treatment of the energy equation at an opening is a combination of that at an inlet and an outlet, depending on whether the flow is in or out.

5.3.5. Symmetry Plane

For a symmetry plane, there is no contribution to the energy conservation equation.

6. Parallelization Strategy

Parallelization is crucial to accelerating the solution process. To describe the communication process between partitions, Figure 6a,b will be used, with different shadings relating to the different partitions. To synchronize the partition, a set of ghost elements is defined at each of the partition’s interfaces. The ghost element adds a layer of nodes to one partition from elements associated with its neighboring partition. This is portrayed in Figure 6c.

It should be noted that the ghost elements are updated in the inter-partition synchronization and are used in order to allow the discretization process to be performed for the nodes that lie exactly at the interface, since they require additional information to reconstruct their missing control volumes. The process starts during the partitioning when the ghost nodes are added to all partition interfaces. During discretization, care is taken to ensure that the matrix coefficients are constructed only for the core elements, i.e., the elements that are not ghost elements. Figure 7 shows the communication process taking place among the different nodes on the partitions. A communication process is, in fact, a synchronization process, through which a value of the node which belongs to its owned partition is copied to its shadow node on the neighboring partition. A shadow node can either be a shared node or a ghosted node; a shared node is a node which is shared by multiple partitions while it does not belong to the new layer of nodes added. Ghosted nodes are simply the nodes that are added to the partition. There will eventually be some useless nodes among the ghosted ones, which are those not connected to any of the stencils; these are also indicated in the figure.

Based on the preceding information, it should be mentioned that synchronization processes are applied at different instants in the solver, some of which are related to fields and others that are related to the algebraic multigrid solver. Flow fields are always synchronized at the beginning of each nonlinear loop so that the assembly process at the nodes associated with the ghost nodes is consistent. For the linear multigrid solver, the parallelization process includes synchronization processes for matrix-vector operations, but there are also multiple requirements that must be satisfied for a fully functioning linear multigrid solver. In this work, the strategy follows the methodology revealed in [36].

7. Test Cases

In this section, the validity of the solver is tested for a set of benchmark problems. Results for each of the following benchmark cases will be provided in comparison with reference data as well as available simulation results from selected commercial packages [37]. Convergence characteristics will also be provided, showing the level of drop of residuals in the flow conservation equations. A representative residual $\tilde{r}$ for a general conservation equation of a quantity $\varphi$ is given by

$$\tilde{r}=\left(\sqrt{\frac{{\sum }_{i,C}{\left(r_i/a_C^{\varphi }\right)}^2}{\overline{n}}}\right)/{\varphi }_{max}$$

(36)

where $r_i$ is the calculated residual value at the node of index i, and is given by

$$\left[r_i\right]=\left[b_i\right]-\left[A_{ij}\right]\left[{\varphi }_i\right]$$

(37)

7.1. Turbulent Flow over a Flat Plate

Turbulent incompressible flow over a flat plate is a standard turbulence validation test case promoted originally by NASA Langley Research Center [38]. Further theoretical and experimental resources of the benchmark are provided by Refs. [39,40]. On this geometry, two meshes were used: a coarse mesh to test the wall function ($35\times 25$) and a fine mesh ($545\times 385$) to investigate low Reynolds modeling. Figure 8 shows a schematic of the flat plate case domain, while Figure 9 shows sample meshes of the case for the two scenarios, the coarse and the fine. Fixed uniform velocity of 69.44 m/s is set at the inlet. Turbulent kinetic energy, as well as turbulence dissipation frequency, are set to 8802.12 (m²/s²) and 126.75 × ${10}^6$ Hz respectively, calculated from a turbulent intensity ratio of 0.039% and an eddy viscosity ratio of 0.009. Convergence characteristics of the two cases are shown in Figure 10a,b, and they appear to be roughly comparable; 150 iterations were sufficient enough to drop the residuals to below 1 × 10⁻⁷. The results of the simulations are stated next. Figure 11a,b show $R{e}_{\theta }$ as a function of x and $C_f$ as a function of $R{e}_{\theta }$ for the two meshes in addition to theoretical data retrieved from [39]. $R{e}_{\theta }$ is defined as follows:

$$R{e}_{\theta }=\frac{{\rho }_{\infty }{v}_{\infty }\theta }{{\mu }_{\infty }}$$

(38)

where $\theta$ is expressed as

$$\theta ={\int }_0^{{\delta }_{99}}\frac{v}{v_{\infty }}\left(1-\frac{v}{v_{\infty }}\right)dy$$

(39)

${\delta }_{99}$ is the boundary layer thickness when the velocity of the boundary v reaches 99% of the free-stream value. $C_f$ is the skin friction coefficient. Distributions of $R{e}_{\theta }$ and $C_f$ are in good consistency with the provided theoretical data. Moreover, Figure 12 shows the velocity wall $u^{+}$ versus non-dimensional distance $y^{+}$ for the two meshes, and results are also consistent with experimental data from Ref. [40].

7.2. Buice 2D Diffuser

The second test case involves the separated incompressible flow through a 2D asymmetric diffuser. The setup is well documented with experimental data available in [41]. The objective of the test case is to prove the ability of the current solver to predict the correct separation behaviors with smooth walls and adverse pressure gradients. A schematic of the case is shown in Figure 13. The parameter H is taken to be 0.015 m. A sample of the corresponding mesh is shown in Figure 14; it includes a total of 133,386 nodes. A velocity of 20 m/s is fixed at its inlet, while a pressure of 0 Pa is fixed at its outlet. Turbulent kinetic energy and dissipation frequency are fixed to 1.5 J/kg and 17,697 Hz, respectively, at the inlet. The case was made to run sufficiently enough to drop the residuals as low as 1 × ${10}^7$. The convergence characteristics are shown in Figure 15. Figure 16 shows $y/H$ versus $10\left|v\right|/v_{\infty }+x/H$ at multiple sections across the diffuser. The origin of the x-positions indicated in the figure is fixed at the start of the expansion section. Figure 17 and Figure 18 show $C_f$ and $C_p$, respectively, versus $x/H$ for the top and bottom walls of the diffuser. Results are in agreement with experimental data. Axial velocity distributions as produced by the commercial package and the current solver at the expansion section are plotted in Figure 19; the distributions show excellent agreement.

7.3. 3D NACA0012 Wing Wind Tunnel Test

The third test case features an incompressible flow past an NACA0012 finite semi-span wing of chord length $c=1.22$ m, fixed in a wind tunnel at ${10}^{\circ }$ angle of attack. Figure 20a,b show the domain of the case and a sample mesh, respectively. The mesh includes a total of 15,114,920 nodes. As for the boundary conditions, the velocity is fixed at 58.81 m/s for the inlet of the wind tunnel, while a pressure of 0 Pa is fixed at the outlet. The boundary which the wing is attached to is set to be a symmetry plane, and all remaining boundaries are no-slip walls. Convergence characteristics for the case are shown in Figure 21. Sample results of the simulation, i.e., longitudinal and cross-flow velocities, are compared against experimental data retrieved from Ref. [42] along two evaluation lines (Figure 22) downstream from the wing at $x/c=0$ and $x/c=0.24$, respectively. Longitudinal velocities at the two evaluation lines are shown in Figure 23a,b, while cross-flow velocities are shown in Figure 24a,b. Both velocities are plotted as a function of the span-wise direction z and appear to be in harmony with the provided data from the reference. A comparison of pressure coefficient ($C_p$) defined as

$$C_p=\frac{p-p_{\infty }}{\frac{1}{2}{\rho }_{\infty }{v}_{\infty }^2}$$

(40)

shows a very good agreement between the commercial package and the current solver (see Figure 25).

7.4. Transonic Turbulent Flow over ONERA Wing

The last test case is the well-known transonic flow over an ONERA wing, used for external flows validation [43]. The mean-aerodynamic-chord length of the wing is $c_{mac}=0.64607$ m and the semi-span length is $b=1.196$ m. The wing geometry and a sample mesh are provided in Figure 26. The mesh includes a total of 250,680 nodes. The freestream mach number, pressure, and temperature are set to 0.8395, 31,500 Pa. and 255.6 k, respectively. Convergence characteristics of the case are shown in Figure 27, and they depict that 50 iterations were sufficient to drop the residuals in all of the flow equations to less than 1 $\times {10}^{-5}$. The resulting pressure coefficient ($C_p$), defined previously in Equation (40), is plotted in the non-dimensional direction $x/c$ at two span sections of the wing, $y/b=0.2$ and $y/b=0.9$, where c is the chord length at the corresponding span position. The span positions are shown in Figure 28. The results are shown in Figure 29 as compared to the commercial software results and experimental data from Ref. [43]. The results are in agreement with the commercial software and show good agreement with the experimental data. Figure 30 shows the pressure distribution over the ONERA wing as produced by the commercial code and the current CVFEM solver.

8. Closing Remarks

In this work, a pressure-based fully-coupled CVFEM solver is presented. The method used departs from the conventional approaches used for the CVFEM, and rather adapts techniques from the CCFVM for resolving the pressure–velocity coupling, and treating the boundary conditions. The discretization of the equations along with the boundary conditions are detailed, and the solver is validated using a set of test cases. The test cases included a turbulent flow over a flat plate which proved the credibility of the wall function implementation as well as the capability for low Reynolds modeling. The two-dimensional diffuser case was set to test the ability to predict the correct separation behaviors with smooth walls and adverse pressure gradients, and the results showed very good agreement with experimental as well as simulation results. The three-dimensional NACA0012 wing case showed consistent longitudinal and cross-flow results with experimental results, and showed very good agreement in pressure coefficient distributions with simulation results. The transonic ONERA case was set to validate the compressibility capturing capability and has shown very comparable results of pressure coefficient profiles at selected section lines with available experimental data as well data from commercial packages. To sum up, the CVFEM method stands today as a strong candidate for being the basis of CFD solvers. This is mainly attributed to the fact that this method combines the essential conservative feature of the FVM method as well as the extremely useful feature of shape functions in FEM, a feature which allows for higher-order solutions without referring to mesh refinements and/or improvements, as is the case for CCFVM. CVFEM today exhibits better documentations as well as more practical implementations, given the availability of useful open-source libraries which are extremely helpful in developing the method.